Since the date today approximates pi (3.14), I thought I would celebrate by deriving an approximation of pi using light sabors.

We start in a futuristic world where urban wildlife have aquired light sabors and taken over Montreal.

To battle these savage beasts, we must learn to make our own light sabor.

set.seed(1.1)
light_sabor_length <- 20
number_of_animal_sabors_in_a_park <- 5000
distance_between_transects <- 40.1

point_1 <- c(-(light_sabor_length/2),(light_sabor_length/2))
point_2 <- c(0,0)
single_sabor <-cbind(point_1, point_2)
rownames(single_sabor)<- c("point_1","point_2")
colnames(single_sabor)<- c("x","y")
single_sabor

##           x y
## point_1 -10 0
## point_2  10 0

Then distribute the army of sabor-wielding critters across the landscape.

First we rotate around the origin to set the instantanious angle of each animal’s individual sabor.

  angle <- runif(1, min=1, max=360)*pi/180
  xy <- as.matrix(single_sabor)
  
  # Rotation
  cos.angle <- cos(angle)
  sin.angle <- sin(angle)
  after_rotation <- xy %*% t(matrix(c(cos.angle, sin.angle, -sin.angle, 
        cos.angle), 2, 2))
  
after_rotation

##              [,1]      [,2]
## point_1  1.100398 -9.939272
## point_2 -1.100398  9.939272

Then we move (a.k.a. translate) the sabor to a location on the landscape where the critter is located. We assume that all the critters are battling vigously and indescriminintly with each other so they are randomly distributed across space and the angles of all the sabors are uniformly distributed.

  # Translation
  translation.x <- runif(1, min=10, max=90)
  translation.y <- runif(1, min=10, max=90)
  translation <- matrix(
    c(translation.x,translation.x,translation.y,translation.y), 2, 2)
  new_coord <- after_rotation + translation
  colnames(new_coord) <- c("x","y")
  rownames(new_coord) <- c("Point_1","Point_2")

Repeat for an entire box of matches

critter_sabor <- function(){
  single_sabor <-cbind(c(-5,5),c(0,0))
  angle <- runif(1, min=1, max=360)*pi/180
  xy <- as.matrix(single_sabor)
  
  # Rotation
  cos.angle <- cos(angle)
  sin.angle <- sin(angle)
  after_rotation <- xy %*% matrix(c(cos.angle, sin.angle, -sin.angle, 
        cos.angle), 2, 2, byrow=TRUE )
  
  # Translation
  translation.x <- runif(1, min=10, max=90)
  translation.y <- runif(1, min=10, max=90)
  translation <- matrix(
    c(translation.x,translation.x,translation.y,translation.y), 2, 2)
  new_coord <- after_rotation + translation
  colnames(new_coord) <- c("x","y")
  rownames(new_coord) <- c("Point_1","Point_2")
  
  return(new_coord)   
}

critter_sabor()

##                x        y
## Point_1 21.93899 84.59096
## Point_2 30.33012 79.15139

for(i in 1:number_of_animal_sabors_in_a_park){
critter_sabor()
}

surface_line <- function(line_height = 0){
  matrix(c(0,100,line_height,line_height), 2, 2, byrow=FALSE)
  }

##      [,1] [,2]
## [1,]    0   20
## [2,]  100   20

one_sabor <- critter_sabor()
one_surface_line <- surface_line(30)
line.line.intersection(one_sabor[1,], one_sabor[2,], one_surface_line[1,], 
                       one_surface_line[2,], interior.only = TRUE)

## [1] 75.73794 30.00000

Make a list of surface lines and a list of matches

list_of_matches <- array(NA,dim = c(2,2,number_of_animal_sabors_in_a_park))
list_of_surface_lines <- array(NA,dim = c(2,2,6))

for(i in 1:number_of_animal_sabors_in_a_park){
  list_of_matches[,,i] <- critter_sabor()
}

line_placement <- seq(0,100, by=20)
for(i in 1:6){
  list_of_surface_lines[,,i] <- surface_line(line_placement[i])
}

crosses <- matrix(NA, 1, 10, byrow=TRUE)

for(i in 1:dim(list_of_matches)[3]){
for(j in 1:dim(list_of_surface_lines)[3]){
  crosses <- rbind(crosses, c(
    line.line.intersection(list_of_matches[1,,i], list_of_matches[2,,i], 
                           list_of_surface_lines[1,,j],     
    list_of_surface_lines[2,,j], interior.only = TRUE), 
    
    list_of_matches[1,,i], list_of_matches[2,,i], list_of_surface_lines[1,,j], 
    list_of_surface_lines[2,,j])
    )
  }
}

crosses <- crosses[which(is.na(crosses[,1]) == FALSE),]
colnames(crosses) <- c("x.cross", "y.cross", "x.p1", "y.p1", "x.p2", 
                    "y.p2", "x.surface.p1", "y.surface.p1", "x.surface.p2", 
                    "y.surface.p2")
head(crosses)

##       x.cross y.cross     x.p1     y.p1     x.p2     y.p2 x.surface.p1
## [1,] 63.08661      20 64.26028 15.53502 61.71804 25.20647            0
## [2,] 14.37017      80 12.61642 83.36074 17.24275 74.49524            0
## [3,] 49.44207      80 44.92515 85.23335 51.45904 77.66313            0
## [4,] 47.22221      20 49.37535 17.30391 43.13500 25.11787            0
## [5,] 24.43356      20 23.31048 15.67796 25.82545 25.35654            0
## [6,] 56.12789      20 53.77911 24.62028 58.31078 15.70602            0
##      y.surface.p1 x.surface.p2 y.surface.p2
## [1,]           20          100           20
## [2,]           80          100           80
## [3,]           80          100           80
## [4,]           20          100           20
## [5,]           20          100           20
## [6,]           20          100           20

pi_estimate <- (2* light_sabor_length * number_of_animal_sabors_in_a_park)/(distance_between_transects * length(crosses[,1]))
options(digits=20)
pi_estimate

## [1] 3.1427417593382642735